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Abstract 

We describe an algorithm for controlling the relative error in the 
numerical evaluation of a bivariate integral, without prior knowledge 
of the magnitude of the integral. In the event that the magnitude 
of the integral is less than unity, absolute error control is preferred. 
The underlying quadrature rule is positive-weight interpolatory and 
composite. Some numerical examples demonstrate the algorithm. 

Keywords: Bivariate integral; Composite quadrature; Interpolatory quadra- 
ture; Cubature; Relative error; Absolute error; Error control 
MSG 2010: 65D30; 65D32; 65G20 

1 Introduction 

We consider the evaluation of 

b u(x) 



I[G{x,y)] = I j G{x,y)dyd. 



a l(x) 

using cubature based on composite interpolatory quadrature, such that 

I[Gix,y)]-Qc [Gix,y)] 



I[Gix,y)] 



(1) 
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where Qc [G (a;, y)] is the composite cubature of G (x, y), e is a user-imposed 
tolerance, and an estimate of I[G{x,yy\ is not known a priori. In other 
words, we seek to control the relative error in the cubature, without prior 
estimation of the integral. The problem is easily understood with reference 
to ([1]): we have 

I J [G (x, y)] - Qc [G (x, y)]\^e\I{G (x, y)] \ , 

so that e\I [G {x^y)]\ is an absolute tolerance. In principle, interpolatory 
methods readily admit absolute error control but, since I[G{x,y)] is not 
known, we cannot impose e \I [G (x,?/)]| as a tolerance. Controlling the rela- 
tive error is appropriate when dealing with integrals of large magnitude; for 
such integrals, absolute error control can be very inefficient. Again, however, 
this presents a problem, since the magnitude of / [G (x, y)] is not known, so 
we do not even know whether absolute or relative error control should be 
applied. 

The algorithm we present here is, in a sense, a 'first-principles' method, 
since it is based entirely on classical concepts relating to interpolatory quadra- 
ture. The list of references [1-9] is our bibliography, and is drawn from the 
estabhshed literature. 

A note regarding terminology: quadrature refers to the numerical approx- 
imation of a univariate integral, and cubature refers to the numerical approx- 
imation of a multivariate integral. Both terms will be used throughout this 
paper. 



2 The Algorithm 

We transform [a, h] to [0, 1] by means of 

x = {h — a)w + a = rriiw + a, (2) 

where x G [a, 6] , w G [0, 1] and mi has been implicitly defined. 
If 

li = min {/ (x) ,u(x)} 

la,b] 

Ui = max {/ (x) , M (x)| 

la,b] 

then the transformation between and [0, 1] is given by 

y = {ui - li) z + h = + li, (3) 
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where y & [Zi, Ui] , z G [0, 1] and m2 has been imphcitly defined. 
As a result of these affine transformations, 



b u{x) 



1 u{w) 




G (x, y) dydx 




G {w, z) mim2dzdw , 



a l(x) 



liw) 



where 



G{w,z) = G {rriiw + a, 'm2Z + h) 
u {rriiw + a) — li 



u (w) 

T{w) 



1712 

I {niiw + a) — li 
m2 



We determine 



M = max < 1, max 

R 



G (w, z) mim2 



(4) 



where R is the domain of integration defined by the transforms ([2]) and ([3]). 
Note that i? C [0, 1] x [0, 1] . Hence, we define 



g [w,z) = 



G {w, z) 17111712 

M ' 



Now, 



U[9{w,z)] -Qc[g{w,z)]\ 

\MI [g {w, z)] - MQc [g {w, z)] \ ^ Me 

G{w,z) 17117712 — Qc G{w,z) 17111712 



^ Me 



G (w, z) 17111712 



Qc G {W , Z) 17111712 



I 


G {w, z) 17111712 


-Qc 


G {w, z) mim2 


I 


G {w, z) mim2 





G {w, z) 17111712 



I [g {w, z)] 



\I[G{x,y)]-Qc [G{x,y)]\ 
I[G{x,y)] 



\I[9 iw,z) 
e 



V\9{w,z)\\ \Qc\g{w,z) 



In the last inequality, we use the fact that changes in variable preserve the 
value of both the integral and the quadrature-based cubature. 
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Clearly, from the last inequality, 



\Qc[9{w.z)]\ 

is an estimated bound on the relative error in Qc [G (x, y)] = MQc [g {w, z)] . 
This estimate is good if Qc [g {w, z)] is accurate which, in turn, is determined 
by the choice of e. 

Now, assuming |/ [G (x, y)] | ^ 1, 

1 M M 



\Qc[g{w,z)]\ \Qc[G{x,y)]\ \I[G{x, 

and, since M is the maximum possible value of |/ [G {x,y)]\ (by construction, 
see dl])), we have 



\Qc[giw,z)]\ 

provided \I [G {x, y)] \ is not substantially smaller than M. For many practical 
situations, this will be the case. However, we have no prior knowledge of 
\Qc [g (w, z)]\ ~ |/ [g {w, z)] \ , so we must be willing to accept the estimate, 
whatever it may be. Obviously, we cannot expect that the relative error will 
satisfy the tolerance e, even if |/ [g {w, z)] — Qc [g {w, z)]\ does. Note that if 

m9{w,z)] - Qc [g{w,z)]\ = e, 

then |Q^[g^^ ^)]| is not merely an upper bound, but is a very good estimate of 
the relative error itself. 

If |J [G (x, ?/)]! < 1, then the relative error could be considerably larger 
than e, particularly if ~ 0, but in this case we favour absolute 

error control (for reasons to be discussed later), and so the relative error is 
not relevant. The quantity Me is an upper bound on the absolute error. 

If the estimate of the absolute or relative error is considered too large, 
say by a factor of rj, then we simply redo the calculation, this time with a 
tolerance of 

£ 

v' 

This refinement is a very important feature of the algorithm, since it enables 
a desired tolerance to be achieved in a controlled manner, even if it requires 
a repetition of the calculation. We are sure that such repetition is a small 
price to pay for a solution of acceptable quality. 
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3 Bivariate composite interpolatory cubature 



Here, we briefly describe bivariate composite interpolatory cubature, includ- 
ing the relevant error analysis. We will consider the effect of roundoff error 
on error control, and offer a criterion for choosing between absolute and rel- 
ative error control. A reasonable degree of familiarity with interpolatory 
quadrature is assumed. 



3.1 The form of bivariate composite interpolatory cu- 
bature 

The composite quadrature that approximates the univariate integral 

G [x] dx 



is given by 

N N 

Qc [G {x)] = c^G (xi) ^hY^^iG {Xi) , 

i=l i=l 

where the Xi are nodes on [a, h] , the coefficients q are appropriate weights, 
/i is a stepsize parameter representing the separation of the nodes, and the 
reduced weights are Wi — Ci/h. 
The bivariate integral 

h u[x) 



J J G {x, y) dydx 



l{x) 



is approximated by 



Qc [G {x, y)]^hY^i[hY ^3'^ (^^' yj^i^ ' 

i=l \ 3=1 ) 

where Vj^i are appropriate reduced weights, yj^i are nodes along the y-axis on 
[Z {xi) , u {xij\ , and ki are stepsizes, with 

M {Xi) - I {Xi) 

Ki — 



2,i 
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Clearly, bivariate cubature is based on univariate quadrature. We can write 
Qc [G {x, y)] = ^Y^ Cj^.G {xi, Vj^i) , (6) 



i=i j=i 



where 



Gj,i — hwikiVj^i- 



3.2 Approximation error 

The approximation error in Qc [G {x)] is bounded by 

A(r) (b-a) /I'maxlG^")! , 

[a,b] ' ' 

where A (r) is a numerical constant particular to the type of quadrature used 
(e.g. Trapezium, Simpson, Gauss-Legendre), and r indicates the so-called 
order of the quadrature. Hence, for bivariate cubature we have 

d'"G 

A (r) (6 — a) max (u (x) — I (x)) ( h^' max ' " ' 



D 



^h^' max 


d'G 


dx''' 





+ (max kl) max 



as an upper bound on the approximation error. The integers A^^i and A^2,i 
in ([H]) can be determined by setting h = max/cj in the above bound, and 
demanding 



^ e 





d'G 




d'G 


^max 


dx^' 


+ max 







A (r) {b-a)D( 



max 



I d'-G I 



+ max 



d'-G 



(7) 



where the various maxima are found on the region of integration. Then 

h — a 



N2A 



h 



U {Xj) - I {Xj 

k 



Furthermore, the stepsizes h and k must be recalculated to be consistent 
with integer values of A^^i and N2^i, as in 

b — a 



h* 
k* 



U {Xi) - I {Xi) 



2A 



(9) 



and it is these stepsizes that are used in Once the stepsizes have been 
determined, the nodes Xi and yj^i can be found. 

This process of computing stepsizes consistent with a tolerance e con- 
stitute absolute error control in bivariate composite interpolatory cubature, 
and is used in the previously described algorithm to find Qc [g {w, z)] such 
that 

\I[9i'w,z)] - Qc [g iw,z] 



€ e. 



It should be noted that our use of max | | + max 



is conservative, 



and could result in smaller stepsizes than is necessary, for the given tolerance. 
However, in these types of numerical calculations it is always better to err on 
the side of caution. Nevertheless, we should be aware that such a conservative 
approach could result in |/ [g {w, z)] — Qc [g {w, z)]\ -C e, so that |Q^[g^^, ^)]| 
overestimates the relative error. Analytically speaking, the approximation 
error is proportional to 



+ 



(e,c) 



where (^, () and (<y9, 0) are points somewhere in the region of integration - 
but since these points are not known, and we cannot be sure of the sign of 



the derivatives, we use max | | + max 



in the error term, instead. 



3.3 Choosing between absolute and relative error con- 
trol 

From ([1]) we have 

\I[Gix,y)]-Qc [Gix,y)]\^E\I[Gix,y)]\, 

so that relative error control is equivalent to absolute error control with an 
effective tolerance e \I [G {x, y)] \ . Replacing £ in ([7]) with e \I [G {x, y)] \ shows 
that, if I J [G {x,y)]\ > 1, h would be larger than if the tolerance was simply 
e, and if |/ [G (a;, ?/)] | < 1, h would be smaller. Consequently, A^"! and A^2,j 
would be smaller or larger, respectively. Smaller values of Ni and N2,i imply 
greater computational efficiency and so, for the sake of efficiency, we choose 
relative error control when |J [G {x, y)]\ > 1, and absolute error control when 
I J [G {x,y)]\ < 1. When |/ [G {x,y)]\ = 1, the two cases are identical. This is 
why we can impose absolute error control on |/ [g {w, z)] — Qc [g {w, z)] \ - by 
our definition of g, I [g {w, z)] is guaranteed to have a magnitude less than 
or equal to one. 
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3.4 Roundoff error 



It is easily shown (see Appendix) that the roundoff error associated with ([5]) 
is bounded by 

4 (6 - a) Dfi, 

where /i is a bound on the magnitude of the machine precision of the finite- 
precision computing device being used, \G (x, ?/)| ^ 1 on the region of integra- 
tion, and the cubature used is based on positive-weight quadrature. In such 
quadrature, all weights are positive; examples of such quadrature include 
the Trapezium rule, Simpson's rule and all types of Gaussian quadrature. If 
the region of integration has unit area, as does I [g {w, z)], then the roundoff 
error simply has the bound 4/i. The roundoff error represents the minimum 
achievable accuracy in the cubature approximation, and is incorporated into 
the error control procedure by replacing the numerator of ([7]) with 

£ - 4 (6 - a) Dfi. 

Clearly, it makes no sense to impose a tolerance smaller than the roundoff 
error. A typical desktop PC has ~ 10~^^. 



4 Numerical examples 

4.1 Example I: relative error control 

We approximate 

I [G {x, y)]= j j e^^Mydx = 1.92660 x 10^ 

1 x2/5 

using Simpson's rule (r = 4, A (r) = . For ease of presentation we show 
all numerical values truncated to five decimals or fewer, although all our cal- 
culations are performed in double precision. The application of the algorithm 
to this example will be described in detail. With the transformations (using 
= 8/5, Zi = 1) 

7z 1 

X = w + l,y = Y + 

> mi = 1,7712 = - 

5 
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the integral becomes 



1 u{w) 



I [G {w, z) 




5 



l{w) 



We find 



u[w\ 



liw) 



M 



7{w + l) 

25 
7{w + l) 

25 



7 
5 
7 
5' 



max 



max 



D 

The stepsize h is given by 



5.07104 X 10^ 
1.67772 X 10^ 

1.57351 X 10^ 

max ( u iw) — I (w] 



18 
26' 



h 



(i^) (1) (i) 
and so, with e = 10^^^. 



max 



+ max 




5.52707 X 10" 



(10) 



A^i = 1810, h* = 5.52486 x 10"^ 

Here, h* is the length of each simpson subinterval (which contains three 
nodes), and there are 1810 such subintervals. Hence, there are 3621 nodes 
Wi on [0, 1] with separation h*/2 (this is the reason for the factor 16 = 2^ in 

m)- 

The stepsizes k* along the z-axis are found from ([S]) and (jH]) for each 
i = 1,2, . . . , 563, and we find 

maxfc* = 5.52706 x 10"^ 

This enables the nodes Zj,i (j = 1, 2, . . . , N2^i) to be found, for each i. As with 
Wi, the spacing between these nodes is k*/2. It must be noted that N2,i could 
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be zero, in which case k* will be NaN {not-a-number in IEEE arithmetic). 
In such cases, it is appropriate to simply set k* =0. 

Composite Simpson quadrature of g {w, z) is performed along the 2;-axis, 
for each i, yielding the 3621 quantities Qc [g {wi, z)] , which have the form 



k* 

Qc [g {wi, z)] = 



g {wi, zi^i) + Ag {wi, Z2,i) + 2g {wi, Z2,i) + Ag {wi, Z4,i) + 
... + 2g (wi, ZN2,i-2,i) + 4^9 {wi, ^jV2,i-i,j) + 9 {wi, ZN^^^i) 



The integer coefficients in this expression are the weights appropriate to 
composite Simspon quadrature. 

Finally, Simpson quadrature is performed over these quantities along the 
w-axis, to give 



Qc [9 {w, z) 



Hence, 



Qc [g (t^i, z)] + AQc [g (^2, z)] + 2Qc [g K, z)] + 
4(3c [9 (w^4, z)] + ... + 2Qc [g {wn,-2, z)] + 
AQc [g (wjvi_i, z)] + Qc [g {wn^.z)] 

3.79922 X 10"^ 
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I[G{x, y)]^MQc [g{w,z) 
The estimate of the relative error is 



1.92660 X 10' 



Qc [9 (w, z) 



= 2.63211 X 10" 



while the actual relative error is 1.47276 x 10 Clearly, the actual error 
is less than the estimate. This is to be expected when using max | | + 



max 



d^'C 



in the computation of h. Obviously, our value for h is conservative 

(^smaller than actually necessary) and so the actual error is smaller than the 
estimate. Nevertheless, as we have stated earlier, it is better to be more 
accurate than necessary and, since the estimate reflects an upper bound, we 
can be sure that the error is no more than 2.63211 x 10~^. If this level of 
accuracy is acceptable, then the result stands. However, if we desire a relative 
error of no more than lO"^*^, say, we simply repeat the algorithm with 



as the new tolerance. This gives 



e 
264 



9.97014 X 10"^^ < 10"^°, 



Qc [9 {w, z)] 

while the actual relative error is 5.35801 x 10~^^ 
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4.2 Example II: absolute error control 



In this second example, the integral 



4 Ix^ 



I [G {x, ; 




sin {xy) 



dydx 



-0.00734 



will again be approximated using Simpson quadrature but, since it has mag- 
nitude less than one, we will see that absolute error control is more efficient 
than relative error control. There is no need for a detailed exposition, as in 
the previous example, and we simply state our results. 
The transformed integral is 

1 u{w) 



I [G {w, z)] 
u{w) 

T{w) 




93 

— sin ((3w + 1) (31^ + 1)) dzdw 
5 



l{w) 



2 {3w + if 



31 



Using M = ^ and £ = 10~* gives 



3w 
31' 

4 



Me 



0.00186 
0.25340. 



0.00136, 



Qc [g {w, z)] 

The upper bound on the relative error is fairly large. The absolute error is 
estimated by Me ; it is clear that, since M is known, e can be chosen so that 

Me equals some desired value. For example, if we seek an absolute error of 
no more than 10~^, we choose e = 5.37 x 10~^, which gives 

Me = 9.9882 x 10"^ < 10"^ 

e 

Qc [g {w, z)] 

with Ni = 1761 (hence, 3523 nodes on the w-axis). Note that achieving this 
tolerance does not require a repetition of the algorithm, since M is known a 
priori. 

On the other hand, to improve the estimate of the relative error to 10"^ 
requires e = 10~^/136 = 7.353 x 10~^, which results in A^i = 2895, and hence, 
more nodes than are needed to achieve the same tolerance in the absolute 
error. This is consistent with our earlier discussion regarding the efficiency- 
based criterion for choosing between absolute and relative error control. 
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5 Conclusion 



We have reported on an algorithm for controlhng the relative error in the 
numerical approximation of a bivariate integral. The numerical method used 
is positive-coefficient composite interpolatory quadrature. The algorithm 
involves transforming and scaling the integral to one that has magnitude 
bounded by unity, and then applying an absolute error control procedure to 
such integral. The relevant scaling factor is then used to find the approxi- 
mate value of the original integral and an estimate of the relative error (if the 
integral has magnitude greater than unity) or absolute error (if the integral 
has magnitude less than or equal to unity) . The calculation can be repeated 
with an appropriate refinement if the estimated error is considered too large. 
The algorithm proceeds in a systematic and controlled manner, and there 
is no need for any prior knowledge of the magnitude of the integral. Two 
examples with Simpson's rule clearly demonstarte the character of the algo- 
rithm. This work extends other work of ours [10], in which we considered the 
control of relative error in the quadrature of a univariate integral. In that 
work, we designated the algorithm CIRQUE, and so we take the liberty here 
of designating the current algorithm CIRQUE2D. 
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A Roundoff bound 

Using ([5j) an ([6]), we write 

Qc[G{x,y)\ = ^Ci{l + fi^^i) 

i=l 

N2,,. 

X ^ fci (1 + Vj^i (1 + fi^j^i) G {Xi, Vj^i) (1 + 

Afi N2,i 

i=i j=i 

where we have indicated the roundoff error in Cj, ki, Vj^t and G (xj, yj^i) explic- 
itly, and we have ignored higher-order terms in the second line. The roundoff 
error T in the cubature is 

Ni N2,i 

^ = YY '^hiG (Xi, Vj^i) + ^^w,^ + l^v,j,i + /^G,i,i) 

i=l j=l 

Ni N2A 

i=l j=l 

where /i is a bound on j| , i \ , \lJ'v,j,i\ \l^G,j,i \ ' ^^'^ have assumed 
\G{xi,yjA)\ ^ 1. Now, since Cj^t = hwikiVj^i = CikiVj^i, 

Ni I N2.i 

T ^ 4/i ^ Ci 1 ^ hvj^i 
i=\ \j=i 
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But, in positive- weight univariate composite interpolatory quadrature, the 
sum of the weights is simply the length of the interval of integration, and so 



T ^ 4/i (6 — a) (max (m (xj) — / (xj))) 
= 4^{b-a)D 



ii{b-a)D^ 1. 
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